Exercise_3.1_Estimating_and_plotting_a_selection_surface_2014.R

josef — Aug 5, 2014, 11:04 PM

## # Computational Exercise 3.1: Estimating and plotting a selection surface
## ### 2012 Stevan J. Arnold

## ## Get your data into R using the following steps.
## Save the file radix5.txt to your desktop or folder of choice
## In the following example, we have created a folder called ?R wd? and placed the data file in it
## Now, set your working directory using a statement comparable to this one
## Be sure and use the / not the \ spacing convention!
#setwd('set working directory here')

## Now read in your data
thamnophis = read.table("radix5.txt",header=T)
#Print out the data frame to make sure it looks OK
head(thamnophis)

  LITTER NAME   SPEED BODY TAIL
1    761    1 0.16207  154   77
2    761    3 0.09019  149   70
3    761    4 0.10823  153   69
4    761    5 0.27609  154   73
5    761    6 0.13532  150   72
6    761    8 0.20325  149   68

## ## Stage 1: Estimating the selection gradients
## Use the following statements to create standardized versions of BODY, TAIL, and TIME
new.body=mean(thamnophis$BODY)-thamnophis$BODY
new.tail=mean(thamnophis$TAIL)-thamnophis$TAIL
new.speed=thamnophis$SPEED/sd(thamnophis$SPEED)
## What do these transformations do to the mean and standard deviation of the variables?
## Find out by using statements like 'mean(new.body), sd(new.body)'.
## Let's begin by fitting a plane to the transformed data
model1 <-lm(new.speed~new.body + new.tail)
## print out the coefficients and statistics of the fit
## this output will give us our best estimates of the directional selection gradients for 
## new.body and new.tail
summary(model1)

Call:
lm(formula = new.speed ~ new.body + new.tail)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.2202 -0.6755  0.0165  0.7030  2.4108 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.62196    0.08393   31.24   <2e-16 ***
new.body    -0.02435    0.02497   -0.98     0.33    
new.tail    -0.00131    0.02800   -0.05     0.96    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1 on 140 degrees of freedom
Multiple R-squared:  0.00684,   Adjusted R-squared:  -0.00734 
F-statistic: 0.482 on 2 and 140 DF,  p-value: 0.618

## Now, fit a full quadratic model to the data, remembering to use the factor of 0.5 for the 
## stabilizing selection gradients
## Create a new variable called prod which is the product of new.body and new.tail
prod = new.body*new.tail
## Then fit the quadratic model
model2 <- lm(new.speed ~ new.body + new.tail + I(0.5*new.body^2) + I(0.5*new.tail^2) + prod )
summary(model2)

Call:
lm(formula = new.speed ~ new.body + new.tail + I(0.5 * new.body^2) + 
    I(0.5 * new.tail^2) + prod)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.2180 -0.6022  0.0418  0.6727  2.4024 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          2.62844    0.11395   23.07   <2e-16 ***
new.body            -0.02759    0.02489   -1.11    0.270    
new.tail             0.00809    0.02800    0.29    0.773    
I(0.5 * new.body^2) -0.00243    0.00914   -0.27    0.791    
I(0.5 * new.tail^2) -0.00171    0.01336   -0.13    0.899    
prod                 0.02031    0.00789    2.57    0.011 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.99 on 137 degrees of freedom
Multiple R-squared:  0.0537,    Adjusted R-squared:  0.0192 
F-statistic: 1.56 on 5 and 137 DF,  p-value: 0.177

## In summary, using z1 for body and z2 for tail,
## we have from the results above, the following values for selection gradients
beta1= coef(model1)[2]
beta2= coef(model1)[3]
gamma11 = coef(model2)[4]
gamma22 = coef(model2)[5]
gamma12 = coef(model2)[6]
## ## Stage 2: Plotting the selection surface, an ISS
## Set the parameters for a full quadratic fitness surface for two traits
z1 = new.body
z2 = new.tail
## Set the value for the intercept of the surface when z1=z2=0
alpha=1
## Define a function called fit that will compute the value of fitness as a function of z1 and z2
fit <-  function(z1,z2, alpha, beta1, beta2, gamma11, gamma22, gamma12)  alpha + (beta1*z1) + (beta2*z2) + (gamma11*0.5*(z1^2)) + (gamma22*0.5*(z2^2)) + (gamma12*(z1 * z2)) 
## Define a series of values for z1 and z2 that will be used to compute the value of relative fitness on the surface
x <- seq(-2, 2, length = 30)
y <- seq(-2, 2, length = 30)

## Compute the surface values of relative fitness using the x-y grid of values for z1 and 
## z2 for later use by the surface plotting function called persp
## We could use the fit function to compute values of fitness but instead we will use a #function, called outer, that is more compatible with persp
## The function outer has some specific requirements so we oblige by writing our fitness 
## function in the following form
z <- outer(x, y, function(a, b, alpha, beta1, beta2, gamma11, gamma22, gamma12) 1 + ( -0.024351*a) + (-0.001307*b) + (-0.002425*0.5*(a^2)) + (-0.001705*0.5*(b^2)) + (0.020315 *(a * b)))

## Define two variables that give the number of rows and cols in z
nrz <- nrow(z)
ncz <- ncol(z)

## Create a function interpolating colors in the range of specified colors
jet.colors <- colorRampPalette( c("yellow", "orange") )

## Generate the desired number of colors from this palette
nbcol <- 100
color <- jet.colors(nbcol)

## Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]

## Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)

## Finally, plot the ISS
par(bg = "white")

persp(x, y, z, col=color[facetcol], phi=30, theta=-30)

## For another view of the surface, do a simple contour plot

par(pty ="s")
contour(x, y, z)

## Let's plot the eigenvectors of the gamma-matrix on this contour plot
## First, let?s define the matrix
gamma = matrix(c(gamma11, gamma12, gamma12, gamma22), nrow=2, ncol=2)
gamma

          [,1]      [,2]
[1,] -0.002425  0.020315
[2,]  0.020315 -0.001705

eigen(gamma)

$values
[1]  0.01825 -0.02238

$vectors
        [,1]    [,2]
[1,] -0.7008 -0.7133
[2,] -0.7133  0.7008

## Next, let?s setup the beta vector and check to see if its values are correct
beta = matrix(c(beta1, beta2), nrow=2, ncol=1)
beta

          [,1]
[1,] -0.024351
[2,] -0.001307

## Using gamma and beta, we can solve for the stationary point on the surface using
## expression(11) from Phillips & Arnold 1989
z.zero = -(solve(gamma)) %*% beta
## Plot the stationary point on the surface
points(z.zero[1,1], z.zero[2,1])
## Finally, plot the eigenvectors on the surface
## First one eigenvector
delx1 = z.zero[1] + 6*eigen(gamma)$vectors[1,2]
delx2 = z.zero[2] + 6*eigen(gamma)$vectors[2,2]
delx11 =  z.zero[1] - 6*eigen(gamma)$vectors[1,2]
delx22 =  z.zero[2]  - 6*eigen(gamma)$vectors[2,2]
x.values=c(delx1, z.zero[1], delx11)
y.values=c(delx2, z.zero[2], delx22)
lines(x.values, y.values, lty=2,lwd=2)
## then the other
delx1 = z.zero[1]  + 6*eigen(gamma)$vectors[1,1]
delx2 = z.zero[2]  + 6*eigen(gamma)$vectors[2,1]
delx11 = z.zero[1]  - 6*eigen(gamma)$vectors[1,1]
delx22 = z.zero[2]  - 6*eigen(gamma)$vectors[2,1]
x.values=c (delx1, z.zero[1], delx11)
y.values=c(delx2, z.zero[2], delx22)
lines(x.values, y.values, lty=2, lwd=2)

## Which of these eigenvectors is gamma max?
## Now, let?s estimate the omega-matrix, approximating it as the negative inverse of the 
## gamma-matrix
omega=-(solve(gamma))
omega

        [,1]    [,2]
[1,]  -4.173 -49.723
[2,] -49.723  -5.936

## It's eigenvectors should be the same as those of the gamma-matrix
## Let's check
eigen(omega)

$values
[1]  44.68 -54.79

$vectors
        [,1]   [,2]
[1,] -0.7133 0.7008
[2,]  0.7008 0.7133

## What is the interpretation of these eigenvalues on the contour plot?
## ## Stage 3: Can you build a loop that will show the contour plot while bootstrapping #over the snake sample? Here?s an example of script that does that!
## We begin by writing a function that we will use during bootstrapping to estimate the coefficients that describe the surface for a particular boot sample
est.coeff=function(y,x1,x2){
  prod=x1*x2
  m1=lm(y~x1+x2)
  m2=lm(y~x1+x2+I(0.5*x1^2)+I(0.5*x2^2)+prod)
  c(m1$coeff[2],m1$coeff[3],m2$coeff[4],m2$coeff[5],m2$coeff[6])
}


## Next, we program a perspective plot function
x <- seq(-2, 2, length = 30)
y <- seq(-2, 2, length = 30)


## Compute the surface values of relative fitness using the x-y grid of values for z1 and 
## z2 for later use by the surface plotting function called persp
## We could use the fit function to compute values of fitness but instead we will use a #function, called outer, that is more compatible with persp
## The function outer has some specific requirements so we oblige by writing our fitness 
## function in the following form
## Write a function that plots the perspective plot for each bootstrap
persp.boot=function(b1,b2,y1,y2,y12,x=seq(-10,10,length=30),y=seq(-10,10,length=30)){
  z <- outer(x, y, function(a, b, alpha, beta1, beta2, gamma11, gamma22, gamma12) 1 + ( b1*a) + (b2*b) + (y1*0.5*(a^2)) + (y2*0.5*(b^2)) + (y12*(a * b)))
  nrz <- nrow(z)
  ncz <- ncol(z)
  jet.colors <- colorRampPalette( c("yellow", "orange") )
  nbcol <- 100
  color <- jet.colors(nbcol)
  zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
  facetcol <- cut(zfacet, nbcol)
  par(bg = "white")
  persp(x, y, z, col=color[facetcol],phi=30,theta=-30)
}

## Test it out
persp.boot(-2,2,1,1,1)

## Bootstrap over individuals
par(ask=FALSE)
n=length(new.body)
for (i in 1:100){
  samp=sample(1:n,n,replace=TRUE)
  boot.speed=new.speed[samp]
  boot.body=new.body[samp]
  boot.tail=new.tail[samp]
  boot.coeff=est.coeff(boot.speed,boot.body,boot.tail)
  persp.boot(boot.coeff[1],boot.coeff[2],boot.coeff[3],boot.coeff[4],boot.coeff[5],x=seq(-10,10,length=10),y=seq(-10,10,length=10))
Sys.sleep(0.1)
}

## ## Stage 4: Let's build a loop that will show the contour plot while bootstrapping over the snake sample while rotating the surface about this vertical axis?
## First, let's make the perspective plot rotate
k=100
for (i in 1:k){
  persp(x, y, z, col=color[facetcol],phi=30,theta=-(30 +(3.3*i)))
}

## To make the plot rotate while bootstrapping we first
## write a new function that plots the perspective plot for each bootstrap
i=1
persp.boot=function(b1,b2,y1,y2,y12,x=seq(-10,10,length=30),y=seq(-10,10,length=30)){
  z <- outer(x, y, function(a, b, alpha, beta1, beta2, gamma11, gamma22, gamma12) 1 + ( b1*a) + (b2*b) + (y1*0.5*(a^2)) + (y2*0.5*(b^2)) + (y12*(a * b)))
  nrz <- nrow(z)
  ncz <- ncol(z)
  jet.colors <- colorRampPalette( c("yellow", "orange") )
  nbcol <- 100
  color <- jet.colors(nbcol)
  zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
  facetcol <- cut(zfacet, nbcol)
  par(bg = "white")
  persp(x, y, z, col=color[facetcol],phi=30,theta=-30+(3.3*i))
}
## Not let's run the bootstrapping during the rotation
par(ask=FALSE)
n = length(new.body)
k=100
for (i in 1:k){
  samp=sample(1:n,n,replace=TRUE)
  boot.speed=new.speed[samp]
  boot.body=new.body[samp]
  boot.tail=new.tail[samp]
  boot.coeff=est.coeff(boot.speed,boot.body,boot.tail)
  persp.boot(boot.coeff[1],boot.coeff[2],boot.coeff[3],boot.coeff[4],boot.coeff[5],x=seq(-10,10,length=10),y=seq(-10,10,length=10))
Sys.sleep(0.2)
}